Harmonic Singular Integrals and Steerable Wavelets 



Here we present a method of constructing steerable wavelet frames 
in L2{M.'^) that generalizes and unifies previous approaches, including 
Simoncelli's pyramid and Riesz wavelets. The motivation for steerable 
wavelets is the need to more accurately account for the orientation of 
data. Such wavelets can be constructed by decomposing an isotropic 
mother wavelet into a finite collection of oriented mother wavelets. 
The key to this construction is that the angular decomposition is an 
isometry, whereby the new collection of wavelets maintains the frame 
bounds of the original one. The general method that we propose 
here is based on partitions of unity involving spherical harmonics. A 
fundamental aspect of this construction is that Fourier multipliers 
composed of spherical harmonics correspond to singular integrals in 
the spatial domain. Such transforms have been studied extensively in 
the field of harmonic analysis, and we take advantage of this wealth of 
knowledge to make the proposed construction practically feasible and 
computationally efficient. 

1 Introduction 

Building upon [27], our purpose in this paper is to provide a systematic 
and practical approach to the construction and implementation of steerable 
wavelets in higher dimensions. The basis of this construction is the theory of 
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Abstract 



singular integral transforms on M"^ with kernels of the form VL{x/ \x\)/ \x\'^^ 
where is a smooth function defined on S"'^^. Properties of such transforms 
were studied by Mikhlin [T3^ , and Calderon and Zygmund [3] . Useful resources 
for this material are the books of Stein and Weiss, I26j . which we shall use 
as primary references. An attractive feature of singular integral transforms 
is their correspondence with Fourier multiplier transforms. For instance, if 
O is a spherical harmonic, then the singular integral transform corresponds 
to a Fourier multiplier that is a multiple of Vt. Furthermore, the spherical 
harmonics (in particular, the zonal spherical harmonics) satisfy symmetry 
properties which make them ideal for applications requiring rotations. 
The two key ingredients for our construction of steerable wavelets are: 

1) an isotropic, band-limited mother wavelet ip that generates a primary 
wavelet frame of L2(M'^); 

2) a finite collection of functions {m„}"™^'' which generate a partition of 
unity on the sphere: 

^ \mn{u)\^ = 1, 

n=l 

where the rUn are purely polar functions; i.e., m„(aj) = mniuj 

The steerable wavelet frame is then generated by the functions 
The limitation on the functions m„ are minimal; however, we shall focus 
on zonal spherical harmonics as they make implementation more amenable. 
Notice that decomposing a signal in this enlarged dictionary of wavelets 
provides more information about the local orientation of the data; mean- 
while, the partition of unity property guarantees that the frame bounds are 
preserved. 

We shall devote the remainder of this section to some basic notation. In 
Section [2j we shall recall some results about singular integral transforms and 
their relation to tight frames. In Section [3j we shall cover some details about 
spherical harmonics and properties of related transforms. Finally, in Section 
|4j we shall describe the steerable wavelet construction, and in Section [5] we 
conclude with some specifics about wavelets based on spherical harmonics. 

1.1 Notation 

The function spaces that we shall consider are the Lp(M'^) spaces, with norm 

II/IImm^) = (£ 1/(^)1" 
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for 1 < p < oo, and our primary focus shall be with L2(M ). Following the 
notation of [SB], we define the Fourier transform of a function / G Li(M"') to 
be 

f{u) = T{f}{uj) = [ f{xy^^dx 
= f /(x)e"2"^^"^da; 

jRd 

and the inverse Fourier transform of / is denoted by J-^^{f}. For a radial 
function f{x) = fr{\x\), we can write the Fourier transform as 



Jo y§d-i \m/ 

/•oo 



{d-2)/2\ 
JO 

where a is the usual surface measure on S"^"^ and J(d-2)/2 is the Bessel 
function of the first kind of order (d — 2)/2, cf. Section VIII. 3]. The 
area of the sphere is given by 



r(d/2)' 

where F denotes the Gamma function. 

For vector valued functions / : M'^ — t- C^, we denote the nth component 
by [f]n- The space of such functions, all of whose components are L2(M'^) 
functions, will be denoted by (K'^), and we define the norm 



N 

L^(Rd) - \ 2^ ll[/]n|lL2 



1/2 



2 Singular integrals and Fourier multipliers 

In this section, we recall some relevant results from the theory of singular 
integrals and provide a basis for the construction of steerable wavelets. One 
of the key ingredients in this construction is a collection of functions which 
generate a partition of unity. In this section we shall show how particular 
classes of such functions behave as Fourier multipliers. 
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Definition 2.1. A collection of complex valued functions M. = {'mn}n=T 
will be called admissible if 

1. Each rUn is Lebesgue measurable and homogeneous of degree 0; i.e., 
mn{aoj) = oPmnioj) = mn{oj) for all a > and u} ^ 0; 

2. The squared moduli of the elements of M form a partition of unity: 

^max 

|mn(c 



n=l 



for every oj G 

The partition of unity property implies that |m„(a;)|^ < 1, so each function 
is a valid Fourier multiplier on L2(M'^); i.e. 



^ II-' 11^2 



Therefore, we can define a transform mapping L2{^'^) to the vector valued 
space i^'^) as follows. 

Definition 2.2. Given an admissible collection M., define the transform 

Tm : L2{R'^) U^"'^-{R'^) by 

[TM{f)]n = J'-\rnnf}, 
and its adjoint T)^ : L^™(M'') ^ L2(M'^) by 

Note that, since we apply these transforms to wavelets, the homogeneity 
condition makes sense, as it means that nin is invariant to scaling. 

Property 2.3. The transform T_\4 maps a wavelet family into another one 
in the sense that 

[TM{H-/a - b))l (x) = [TMmn - b) 
for any ip G L2(M''), a G M+, and b G R'^. 



Also, it follows from Plancherel's identity and the partition of unity 
condition that Tv; is in fact an isometry. Hence, we can apply this transform 
to a tight frame to generate a new frame with the same frame bounds. 
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Theorem 2.4. Suppose {(j)k : A: G Z} is a Parseval frame for L2(M'^); i.e., 



and 



W 

k 



for every f G L2{W^). Then the multiplier transform Tvj associated with an 
admissible collection generates a related Parseval frame: 

{'0fc,n = [TM{(t>k)]n : G Z,n = 1, . . . ,nmax} , 

with 



n=l k 

Proof. Our proof follows the same lines as [29i. Proposition 1]. Notice that if 
is admissible, then so is M := {mn}"™. From the definition of Ai 

Additionally, each component of Tj^f can be expanded in the original frame 

•' M 

where F is the function with components 

[F]n = Y.(\-^M(f)]n,cl>k)<Pk 
k 

= "^{f^ [TMiMU) (pk 



The reproduction property now follows by computing the product T^F: 

^ = -^^ 1 E E i'^Mmu) ^k \ 

[ n k ) 

n k 
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To verify that the frame is stih tight, write 



\iTM{f)]n\\; = J2\{[Tj^{f)]n, 
k 



k 

SO that summing over n gives the result. □ 



Notice that our Definition 2.1 of admissibility is fairly general and directly 
exploitable for implementing wavelets in the Fourier domain. However, if a 
feasible spatial domain representation is required, we must impose certain 
restrictions. A reasonable condition is to assume that the elements of M are 
smooth, since it allows us to relate the multiplier transform to a singular 
integral transform. 

Theorem 2.5. 124[ Theorem III. 6] Let m be homogeneous of degree and 
indefinitely dijjerentiable on S'^^^. Then for 1 < p < oo the Fourier multiplier 
transform T : Lp{M.'^) — )■ Lp{M.'^) given by 

T{f) = F-\mf} 
can be computed by the singular integral 

T{f){x) = cfix) + hm / ^fix - y)dy, 

where c is the mean value of m on S'^^'^ C M'^ and the functions m and Q 
are related by 



m{u) = c + j^^ ^ sign(cj • y) + log j^^^T^ J J ^{y)<^CF{y). 

3 Spherical harmonics and singular integral oper- 
ators 

The class of smooth admissible functions with which we shall be primarily 
concerned are the spherical harmonics. These can be viewed as a multi- 
dimensional extension of the trigonometric polynomials in one dimension. 
The key property is that it is possible to construct an orthonormal basis of 
L2(S'^~^) using spherical harmonics. As Fourier multipliers, the spherical 
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harmonics of degree one correspond to the Riesz transform, which has 
previously been used to construct steerable wavelets It turns out that 
transforms based on any of the spherical harmonics have similar properties. 
One remarkable aspect is the symmetric role assumed by the spatial and 
frequency variables and the fact that the Fourier transforms are part of the 
same family, cf . [26', Chapter IV] , [24', Chapter III] . For the benefit of the 
reader, we review the properties of the spherical harmonics that are relevant 
for our purpose. 

To begin, we recall that a homogeneous polynomial of degree i on M"^ is 
a linear combination of monomials 

\ol\=1 

where the degree |q| = ai + • • • + of each monomial cj" = o;"^ • • • uf^^ is 
i and each Ca is complex. A homogeneous polynomial P is harmonic if it 
satisfies Laplace's equation 




The spherical harmonics of degree ^, denoted by can then be defined 
as the restrictions to the sphere of such polynomials. Specifically, every 
homogeneous harmonic polynomial P of degree I defines a spherical harmonic 
Y S M'n by the equation 

Y{u) = P{u,) 

for any u G C M'^. Moreover, any homogeneous polynomial P of degree 
/c > can be expanded as 

P(a;) =Po(^) + --- + |^|''i^£H, 

where each Pj is a homogeneous harmonic polynomial of degree k — 2j, 
cf. |26| Theorem IV. 2.1]. This means that the proposed steerable wavelet 
construction includes previous constructions that utilized higher-order Riesz 
transforms. 

An important property that is implicit in many results is that spherical 
harmonics of different degrees are orthogonal on the sphere S*^"^, cf. \14:\ 
Lemma 2] or [26. Corollary IV. 2. 4]. In particular, since Pq is constant, this 
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implies the spherical harmonics of positive degree have mean zero on the 
sphere: 

/ Y{uj)da{u) = 

for any Y e J^i with ^ > 0. 

The dimension of Jifg can be computed to be 

We shall use the following notation to denote a real-valued orthonormal basis 

\Ye,k:S''-^ ^C:k = l,...,Nid,£)]. 

Summing over i < imax, we can determine the dimension of the space of 
spherical harmonics of degree at most ^max to be N{d + l,^max), cf. [U p. 
4] or [31, Chapter 17]. Explicit constructions of these basis functions are 
known. For example, the three-dimensional case is analyzed in detail in [21 
Chapter 9], and a general approach for higher dimensions is given in 
Chapter 2]. 

We conclude this review with the property that is most important for 
our purpose: any orthonormal basis {Yi^k} of J^i satisfies 

cf. \14:\ Theorem 2] or |26| Corollary IV. 2. 9]. This is indeed a powerful 
property, for it implies that the collection of Fourier multipliers 

is admissible. Furthermore, as we shall see shortly, this property allows us to 
define a collection of admissible multipliers which contains all of the spherical 
harmonics up to a fixed degree £max- 

3.1 Harmonic Riesz transforms 

The motivation for steerable wavelets is to provide a wavelet decomposition 
which more accurately accounts for local orientation of data. In a series 
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of papers [271 EH [22] , steerable wavelets have been constructed using the 
Riesz transform and its higher-order variants, taking advantage of their scale 
and rotation invariance and their unitary character. We recall that the 
Fourier multiplier of the order i Riesz transform of an L2iW^) function / 
is a vector valued function whose components are (-y/£!/Q;!)cj"/ where 
|q!| = ^. In the spatial domain, this translates into a principal value singular 
integral. The use of spherical harmonics generalizes but also simplifies this 
construction, thanks to its orthogonality properties. 

To make our construction precise, we designate the Fourier multiplier 
transforms associated with the spherical harmonics as harmonic Riesz trans- 
forms. 

Definition 3.1. For any positive integer ^max and any unit vector c = 
(co, . . . , ci^^^) G M^™=''=+^, we define the order 

^max harmonic Riesz transform 

to be the multiplier transform Tvj, where 

where i ranges from to £max and k ranges from 1 to N{d,i). 

Note that if c contains entries which are zero, the corresponding multi- 
pliers are not included in the transform. The admissibility of this transform 
follows immediately from Equation ([3]). Furthermore, when applied to smooth 
functions with vanishing moments, the transform preserves decay as well as 
vanishing moments. 

Theorem 3.2. Let ijj he a differentiable function (or wavelet) with vanishing 
moments of order N > 1 such that tp and its derivatives satisfy the decay 
estimates 

1. \^j{x)\ <C{l + \x\)-'^-^+', 

2. \D°'ip{x)\ <C{l + \x\)-'^-^-'^+^, \a\ = l 

for some C > and < e < 1. Then for any i >0 and any 1 < k < N{d,i), 
the corresponding component of any harmonic Riesz transform Tj^('iIj) has 
decay similar to ^ and maintains the same number of vanishing moments, 
i.e. 

\[TMmeM^)\<C{l + \x\r^-^+^' 
for some < e' < 1 and [Tj^{ip)]i^k has N vanishing moments. 
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Proof. This follows from the proofs of |30| Theorems 3.2 and 3.4]. Those 
results were stated for the first-order Riesz transform; however, one can verify 
that they hold for a more general class of principal value singular integral 
operators. 

The kernels of the singular integral operators that define the Riesz 
transform are K{x) = x^^j ^ while the kernels used in the harmonic 

Riesz transforms have the form K = P{x)/ \x\'^^^ , where P is a homogeneous 
harmonic polynomial of degree I. The essential properties of Riesz kernels 
that were used in the proof of those theorems were: 

• K has mean zero on the unit sphere; 

• the kernel K and its derivatives satisfy certain decay conditions; 

• K \s smooth away from the origin, so that it can be well approximated 
locally by polynomials. 

It can be verified that each of these conditions holds for and hence the 
results of [30] are applicable as well. □ 

3.2 Generalized harmonic Riesz transforms 

While directly using the spherical harmonics in a steerable wavelet frame 
provides a means of categorizing data, we would like to extend this method 
to make our wavelets adaptable and possibly easier to implement. The 
approach we take is to compose the harmonic Riesz transforms with matrices 
representing isometrics. Using an isometry, the derived collection of Fourier 
multipliers will again be admissible. 

Definition 3.3. Let T_m be a harmonic Riesz transform consisting of N 
elements. Additionally, let U be a complex valued matrix of size rimax x 
which represents an isometry; i.e., U^U is the identity matrix of size A^. We 
define the associated generalized harmonic Riesz transform of / G L2(M'^) 
to be the vector- valued function Ta4,u(/) obtained by applying U to the 
harmonic Riesz transform Ta/((/) of /; i.e., the components of Tv(,u(/) are 
linear combinations of the components of Tj^{f). 

Note that in this definition, we deal with two admissible families. There- 
fore, we have denoted the size of the original family as N, in order to reserve 
'^max for the size of the derived family. Also, note that a consequence of the 
isometry condition is that rimax ^ A^- Additionally, if the initial harmonic 
Riesz transform Ta/j is defined by a vector c with all non-zero entries, then 
in the above definition A^ = N{d + l,^max)- 
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4 Steerable wavelets 



In previous sections we covered the mathematical tools necessary to transform 
a wavelet frame into a steerable one. In this section, we complete the 
construction by introducing an appropriate primal wavelet basis. Our choice 
is the direct extension of the two dimensional case [27j . 

Proposition 4.1. Let h : [0, oo) — )• M 6e a smooth function satisfying: 

(1) h{oj) = for > 1/2 

(2) 5]|M2M|' = 1 



(3) ■l""*'-' 



do;*^ 



= Oforn = 0,...,N. 

uj=0 



Then the isotropic mother wavelet ^ whose d-dimensional Fourier transform 
is given by 

4^{u) = h{\u\) 

generates a tight wavelet frame o/L2(M°') whose basis functions 

tpj,k{x) = il^jix - 2^k) with i/jjix) = 2~^''^/V(2"-''«) 

are isotropic with vanishing moments up to order N . Additionally, any 
L2(M'^) function f can be represented as 

Proof. This follows from a combination of Parseval's identity for Fourier 
transforms and Plancherel's identity for Fourier series. □ 

We now define a steerable wavelet frame to be a generalized harmonic 
Riesz transform of a primal isotropic frame. Theorem |2 .4| guarantees that the 
frame bounds are preserved and that we maintain the reproduction property: 

j& k&'Ld n=l 

As we are applying the harmonic Riesz transforms to isotropic functions, 
the transform can be reduced to a more manageable form. Essentially, the 



following is a simplification of Theorem 2.5 which uses radial symmetry to 



reduce the Fourier transform to a one dimensional integral. 
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Theorem 4.2. JM Theorem IV. 3. 10] Suppose d > 2 and i/j e L'^{R'^) n 
L-^(M'^) has the form 



where P is a homogeneous harmonic polynomial of degree i, then ip has the 
form ip{x) = F{\x\)P{x) where 



Jo 

and Jy is the Bessel function of the first kind of order u. 

5 Directional wavelets using zonal harmonics 

We use the term steerable to convey the fact that the wavelets we construct 
are intended to be rotated (or steered) to provide a better analysis of the 
data. To see how this is accomplished, let us consider a generic steerable 
wavelet ipGen] i-e., an element of the generalized harmonic Riesz transform 
Tm,u(V') where tp is a, primal isotropic wavelet. Each such function is of the 
form 



1=0 k=l ^' 
where the coefficients U£^k are related to the rows of U. Of particular interest 
are the cases where 



for some G ]R'^\{0}. The resulting spherical function is a zonal function, 
and it has the form 



where Pe{d; •) is a generalized Legendre polynomial of degree i, cf. 
Section 1.2]. A formula for these polynomials is 



^P{u) = h{\u\) 



Piu) 
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for X G [—1,1]. 

One benefit of using these zonal functions is that we can compute rotations 
fairly effortlessly using the formula: 

Pi,{d; Ra;o • Rt<j) = Pi{d; uq-ud), 

for any matrix R satisfying R^^ = R-^. In addition to making rotations 
straightforward, this structure implies that the value of V'Gon is determined 
solely by the distance of o^o/ I'^ol from tij/ on the sphere; i.e., it is a zonal 
function. 

Zonal functions have proved to be particularly valuable for approximation 
on spheres. For example, positive linear combinations of generalized Legendre 
polynomials are positive semi-definite functions, which can be used for 
interpolation |THl [Ml EH E2] • Indeed, the steerable wavelet construction we 
propose uses polynomials which are positive semi-definite on the sphere. 



5.1 Two dimensions 

We would now like to present some information regarding the implementation 
of steerable wavelets based on generalized Legendre polynomials, and we shall 
start with a two-dimensional motivating example. A basis for the circular 
harmonics of degree I on is 

{sin(£a;), cos(^(^)}. 



Considering Definition 3.3, we shall define an isometry to generate a new 
partition of unity. Note that we shall neglect certain scaling factors, as they 
do not impact the underlying principle. 

Given a point G S^, we define the kernel P in the span of the degree i 
basis by 

P(a;o, w) = sin(^a;o) sin(^a;) + cos{lujo) cos(£u;) 
= cos(^(a;o — uj)). 

For a collection of points {wn}^™^" the matrix that transforms the degree £ 
basis {sin(^ti;), cos(^ti;)} into {cos(£(a; — wi)), . . . ,cos(^(cj — t^^nmax))} is 



/ sin(^a;i) 
sin(^a;2) 



cos{£uj2) 



\sin(^a;„^^J cos(£a;„^^J/ 
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To ensure that the columns of are orthogonal, we require 



= > sm{£ujn) cos{£uJn) 



n=l 

'^max 



= sin(2£a;„). 



n=l 



We could consider choosing the points cOn to be roots of sin(2£-); however, 
this approach would be less useful in higher dimensions. Instead, we shall 
use a circular quadrature rule with equal weights. Let {uJn}n=T ^ 
points for which 



would imply that the columns of XJi are orthogonal, and hence that is an 
isometry (after normalization). 

Sets of points satisfying Equation Q are well known and are referred 
to as spherical t-designs, where t indicates the maximum degree polynomial 
for which quadrature holds. Such sets are known to exist for arbitrarily 
large t [21], and the most natural choices consist of equidistributed points 
[1]. Incidentally, Simoncelli's two dimensional equiangular steerable wavelet 
construction can be reinterpreted in terms of such i-designs jl8| . 

5.2 Higher dimensions 

Fix i > and let {Yi,k}^}^'^^ be an orthonormal basis for the spherical 
harmonics of degree £ on S*^"^. As in the two-dimensional case, we shall 
define an isometry to generate a new partition of unity. Specifically, we select 
\Ji to be the matrix satisfying 

for a collection of points X = {un}n=T S"^"^ to be specified. Delsarte 
et. al refer to this as the ith. characteristic matrix associated with X [H 





for all trigonometric polynomials of degree at most 2i. Then 



= 
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Definition 3.4]. Applying to the basis generates a new basis of zonal 
polynomials 



^ lf,fe(a;m)Y£,fc(cj) = -—^-^-Pi{d;u;m ■ uj)- 
k=l ' 

In order to make the columns of orthogonal, we need 

^max 

= ^ Y£,fc(^„)Y£,fc/(cj„) (5) 

n=l 

for k ^ k' . As in the two-dimensional case, we choose the points X to form 
a spherical 2£-design, so that 



Tl 



V p{un) = / p{uj)da{u) (6) 



for spherical harmonics of degree at most 2£. In fact, the conditions on X 
considered in ([s]) and ^ are almost equivalent, cf. [H Remark 5.4]. On S^, 
examples of t-designs are provided by the vertices of platonic solids: the 
vertices of an icosahedron or a dodecahedron constitute 4-designs [8j. More 
generally, t-designs appear as orbits of elements of S"^"^ under the action of 
a finite subgroup of the orthogonal group on the sphere [1]. Some specific 
examples are given in an online library of t-designs f9]. In particular, this 
library contains t-designs on S'^, where t ranges from to 21. 



5.3 Localized kernels 

Concerning the construction of zonal basis functions, one final point to 
address is localization. The reason for this is that well localized functions can 
be used to more accurately represent the orientation of data. While we would 
ideally like to use locally supported functions, they cannot be represented 
as polynomials. Therefore, we shall instead use a normalized polynomial 
approximation of the identity. 

Let us first recall that on the circle, the delta function can be represented 

by 

1 I °° 

'^o(w) = — + - V'cos(^a;), 
zvr TT ^ — ' 

and we could construction an approximation by truncating this series, pro- 
ducing a Dirichlet kernel. The problem with such a construction is that 
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the Dirichlet kernel is highly oscillatory. Therefore, we instead propose to 
construct approximate identities analogous to the scaling functions of Freeden 
et al. [H Section 11.1.3]. For this construction, we begin with a compactly 
supported function a : [0, oo) — )• [0, 1] satisfying: 

1. a(0) = 1, a(l) = 0, and a{uj) > for a; G (0, 1); 

2. a is monotonically decreasing; 

3. a is continuous at and piecewise continuous on [0, 1]. 
Then an approximation to is given by 



costw 



Now, based on the results of the previous section, we know that an 
admissible zonal basis is given by 



k{uj - uJn) = Co J h ^ ceJ cos(^(a; - Un)) 



for any unit vector c and any circular 2^max-design {ujn}n=T- Therefore we 
should choose c to be a normalization of the coefficients of Si^^^ to obtain a 
well localized basis. 

Several examples of such functions appear in the literature. For instance, 
in [6l Section 11.4.3] the authors propose the cubic polynomial 

= (1 -w)2(i + 2w) (7) 

to produce a function Si^^^{uj) with suppressed oscillations. An alternative 
choice, based on the localization analysis of |T6| and a construction from [llj, 
is to choose a to be a B-spline centered at 0. B-splines also appear implicitly 
in Simoncelli's steerable pyramid construction [18j . The angular part of his 
wavelets are of the form 



■"-max 

cos(l^)^— = e'^'-- Yl ^ '"'Me-*"2^ 



4 



i=0 



E 



/ ^max \ „ 

e 



and the binomial coefficients are the discrete analog of the B-splines. 
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In higher dimensions we can apply the same analysis to construct localized 
kernels. This leads us to approximate the delta function at ojq by 




Additionally, our analysis requires that we use kernels of the form 



^ max 




as before, we can adjust the coefficients of Si^^^^^ to determine an appropriate 
c in A. Note that the choice of q and 'n-j^^x 

must be balanced to produce a 
good kernel. Choosing a smoother a produces a kernel with less oscillation; 
however, its main lobe will be less localized. To compensate, we could 
increase the degree £max) but this means that we need to increase nma.x', i-e., 
utilize a larger collection of basis functions. Example zonal kernels A are 
plotted in Figure [T] 

If localization is of primary importance, one can define a variance on the 
sphere, which should be minimized by a polynomial of a given degree. Indeed, 
this approach was used in |27j for steerable wavelets on M?. Furthermore, 
an uncertainty principle for localization in both space and frequency was 
studied in [6l[7l[I7j. Related results are contained in [5j, where the author 
considers polynomials whose degrees lie within a given range. Additional 
work concerning the localization of spherical Slepian functions can be found 



In the two-dimensional case, Unser et al. extend the definition of variance 
and allow for a general class of weight functions [271 Appendix A] . A similar 
generalized definition of variance was introduced by Michel on ^ 
However, instead of using Fourier analysis and Bochner's theorem, his results 
utilize the theory of orthogonal polynomials. In Appendix A, we provide an 
extension of this type of approach to the higher dimensional setting that can 
serve as a basis for the design of wavelets with an optimal angular selectivity. 

6 Construction and steering 

In this section we lay out the construction of the zonal basis and spherical 
harmonic basis in more detail. Furthermore, we show how the wavelets can 
be steered using matrix multiplication. To begin the construction, we choose 
a maximum degree ^max and a real orthonormal basis of spherical harmonics 



in [221123]. 
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Figure 1: Plot of zonal kernels A(cos(-)) for d = 3, imax = 10, and rimax = 216. 
The kernel on the left was constructed using a from Q , while the functions 
on the right were constructed by choosing a to be B-splines. In the plot on 
the right, the dashed line corresponds to a linear B-spline, and the solid line 
corresponds to a cubic B-spline. 
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of degree at most £max in a vector [Y]m = Ym of length N{d + 1, iraax)- Then, 
given a unit vector c = (cq, . . . , Cmax) G M^™^"'^^, we define a diagonal matrix 
C of size N{d + 1, ^max) X N{d + 1, ^max) as follows: For [Y]m of degree i, 
we set 



N{d, i) 



This implies that the entries of CY form an admissible collection. Next we 
choose a collection of points on the sphere Xq = {uJn}n=T^ which form a 
2^max-design on S*^"^. To construct the final admissible collection, we define 
the nmax X N{d + 1, ^max) matrix Uxo as 

[Uxo]n,m = \ ^m(<^n)- 

y '^max 

Our admissible collection is then given by := Ux„CY. 

For steering, we use the isometry property of Uxq • Let R be a rotation 
matrix and define Xi to the collection of points obtained by rotating the 
elements of Xq by R. We then can expand the elements of the rotated basis 
ZiXi in terms of the original basis Zxq as 



Zx,=Vx,CY 

Since the steerable wavelet frame is a Parseval frame, the steering matrix, 
which transforms the wavelet coefficients from the original basis into the 
coefficients of the rotated basis, corresponds to the matrix mapping Zxq 
to Zxi- Hence the rimax x ''^max steering matrix is given by S := UxiU^^^. 
Using the properties of zonal harmonics, we can see that the entries of this 
matrix are pointwise evaluations of a zonal polynomial A^^,^^ (see Figure [2] 
for an example). Precisely 
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Figure 2: Plot of the kernel A£j^^^^(cos(-)) for ^max — 10 and Hmax — 

216. 



_^ Ar{a!+l,^niax) 



m=l 

d—l\ *^max 



^max 



i=0 
and 



Interestingly, the steering operation is very much akin to an interpolation 
that uses A/>„„ kernel (cf. Figure^. 

Depending on the number of points, the steering matrix can be quite large. 
As an alternative, one could work directly with the orthonormal spherical 
harmonic basis. In this situation, the steering matrix reduces significantly as 
it is a block diagonal matrix with blocks of size N(d,£). To see this, let us 
define the vector Y by 
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/ Yo \ 
Yi 



\Y£ / 

where each is a vector whose components are an orthonormal basis of 
the spherical harmonics of degree t: 

In order to make the entries of Y an admissible collection, we multiply by 
the block diagonal matrix 



/Co 
Ci 



\ 




VO ••• 

where is the N{d, I) x N{d, t) diagonal matrix with entries 



(4iax + l)iV(d,^)' 



With C defined in this way, we are giving equal weight to each degree L 
One could choose an alternative weighting, but it must be constant over a 
given degree for the partition of unity property to hold. Now to construct a 
steering matrix, we need to find an expression for the spherical harmonics in 
terms of any rotation of them. Orthogonality between degrees of spherical 
harmonics makes this especially nice because it means that we can rotate any 
given degree independently of the others. Therefore, let us fix i and consider 
the problem of steering the functions in C^Y^. As is a constant multiple 
of the identity matrix, this is equivalent to finding a steering matrix for Y^. 
Since we are dealing with an orthonormal basis, we can expand any as 



N{d,l) 

fe=i 
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and computation of these inner products can be handled using a quadrature 
rule. Specifically, for any 2^-design {i0n}n=ij we have 

n=l 

Therefore, the N{d,£) x N{d,i) matrix with entries 

[^i]k,k' = ^l£,fc'(Ru;„)y£,fc(a;„) 

n=l 

transforms the original basis Yg into the rotated basis Y|^, i.e. V^Y^ = Y|^. 
Consequently, the steering matrix which transforms the wavelet coefficients 
corresponding to Y^ into the coefficients corresponding to Y|^ is also given 
by V^. Note that in two dimensions N{d,i) = 2 for i > 1, and in three 
dimensions N(d,i) = 2i + 1. Therefore, these steering matrices can be 
significantly smaller than the ones used in the zonal construction when a 
large £ma.x is chosen. 



7 Conclusion and practical summary 

Throughout the course of this paper, we have developed the theory of steerable 
wavelets in any number of dimensions greater than one. The previous two 
sections were devoted to several technical aspects of implementation, and 
here we summarize the construction. Since the zonal construction is perhaps 
more tractable conceptually, we now concentrate on this. One can think of 
this construction as a generalization for d > 2 of Simoncelli's equiangular 
design. 

The substructure of a steerable wavelet frame is an isotropic mother 



wavelet ip, which satisfies the conditions of Proposition 4.1 A new, expanded 
frame is produced by the collection of mother wavelets 



where M is an admissible class of functions. Now given a maximum degree 
^max, a unit vector c = (cq, . . . , Q^ax) ^ M^™'"'+-^, and a spherical 2£inax-design 
X = {iUn}^'!^^ define the admissible collection 



M 



mn{u) = } cn Pe, d; r • —r 



: n = L Ur, 
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where the Pe{d; •) are generaUzed Legendre polynomials. The main properties 
of this basis are summarized in the following theorem. 

Theorem 7.1. The wavelet construction given above defines a tight curvelet- 
like wavelet frame of L2iW^), where all wavelets are rotated versions of a single 
wavelet template (per scale). Furthermore, these wavelets are parametrized 
with a set of coefficients q that can he chosen arbitrarily. 

Ideally, the vector c is chosen so that the functions rrin are well localized 
and peaked at u; = u^n- For large collections of points X = {un}^^i on 
the sphere and well-localized functions runii^), this new frame can detect 
the orientation of data and provide more information concerning structure. 
However, one can also work with a smaller collection Ai and take advantage 
of the steering property to orient the wavelet frame in a data-adaptive fashion. 
In order to steer the wavelet basis, we need a matrix which transforms the 
wavelet coefficients upon rotation of the wavelet basis. The structure of the 
zonal basis allows us to use the nmax x '?^max steering matrix 

for a rotation R, where 

1=0 "'"'^^ 

Notice that the entire construction is in terms of the zonal basis, which 
means we can avoid working directly with the spherical harmonics. On the 
other hand, working with the spherical harmonics would allow for smaller 
steering matrices. 

A Optimal spherical polynomial constructions 

Using the work of Michel [T^] as a starting point, in this appendix we 
shall show how orthogonal polynomials can be used to construct spherical 
polynomials that minimize energy functionals on S*^"^. In order to simplify 
computations let us fix a point cjq G S'^~^ and assume that each of the 
generalized Legendre polynomials in {P£(d;uJo ■ <^)}£^'^ has been normalized 
to have -L2(S'^~^) norm one. 

Considering the framework of our problem, we shall address the problem 
of finding the polynomial (centered at uq) 
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£=0 



of norm one in L2 (S ) that minimizes (or maximizes) an energy functional 
of the form 

E{P,;W)= [ \P,{uo-u)\''W{uo-u)daiu), 

where W is an arbitrary positive continuous function. One possibihty would 
be to choose W{t) = arccos(f)^, so that the energy functional would be 

E{Pc, arccos(-)^) = / \Pc{uiQ ■ u>)\'^ arccos(ci;o " <^)'^da{u) 
= / |Pc((^o • ^)|^dist(a;o,t<^)^d(T(a;), 



where dist refers to the spherical distance. 

As the energy functional contains only zonal functions, it can be reduced 
to a simpler form: 



v^r(— ) Jo 

\p^{t)fw{t){i-t'^Y'^-^y^dt. 



2 



-1 



Now notice that any polynomial of degree at most ^max is a linear combination 
of generalized Legendre polynomials {-P^(d; •) : £ = . . . ,£max}- Also, the 
assumptions on W imply that there exists a sequence of polynomials {QiYt^o 
that are orthonormal with respect to E, 



2 



1 



where each Qi is a polynomial of degree cf. [10]. Hence, there exists an in- 
vertible change of basis matrix from {Pi{d] ■)Yi^Q to {QiYi^q • Consequently, 
the polynomial Pc that minimizes the functional E[P(.; W) is determined by 
setting c to be the unit eigenvector corresponding to the minimal eigenvalue 
(in absolute value) of the change of basis matrix. 
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